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Abstract 

The classical SU(3) gauge theory is shown to be deterministic chaotic. Its largest 
Lyapunov exponent is dertermined, from which a short time scale of thermalization 
of a pure gluon system is estimated. The connection to gluon damping rate is 
discussed. 



In Ref. [1], a classical lattice gauge method was introduced showing that SU(2) gauge 
theory is chaotic in the semi-classical limit. The largest Lyapunov exponent was obtained. 
In this Letter, we use the same method to study SU(3) gauge theory. While the study 
of SU(2) theory gives some hints on properties of a general non-abelian gauge theory, 
the present study of SU(3) will give us some direct and specific information about a 
highly excited gluon system. The reason why such a semi-classical approach is relevant 
here is that in a highly excited system the contribution of quantum fluctuations to some 
observables (for example, Debye screening mass, gluon damping rate, we will explain this 
briefly in the end of this Letter) is negligible compared to those from thermal excitations. 
So at least for some observables, classical calculations can give correct answers. Here 
we want to use this method to study entropy production and thermalization in a highly 
excited gluon system. This process is relevant in the study of the early universe 
and of relativistic heavy ion collisons ^ §], where the phase transition of nuclear 
matter to quark-gluon plasma is expected to occur. The proposed signatures of this 
transition depend critically on the time scale of the pre-equilibrium thermalization process. 
Studies ^ show that the time scale of this process is short and the gluon sector 
thermalizes much faster compared to the quark sector (r^ = Sr^). This later fact makes 
it plausible to consider the thermalization of the gluon sub-system separately. But up 
to now, this process is still not completely understood because a fully non-perturbative 
quantum mechnical calculation is beyond our ability. On the other hand, our method 
not only leads to an exact and non-perturbative solution to this problem (in the classical 
limit), but also helps to explain the reason of a fast thermalization process in the classical 
language. That is the chaotic nature of the non-abelian gauge theory. We hope that our 
results will be helpful to a full quantum treatment in the future. 

This Letter is organized as follows. First we will outline the method of a classical 
Hamiltonian lattice simulation. Some technical points will be explained. Then we will 
proceed to the results. An exponentially increasing distance between two field config- 
urations is observed. From this behavior, an energy dependent Lyapunov exponent is 
obtained which leads to an estimate of thermalization time scale for the gluon system. 
Finally we will discuss the relation between Lyapunov exponent and gluon damping rate. 

The framework of a classcal lattice simulation in Hamiltonian formalism of SU(2) 
gauge fields can be found in Ref[l]. Here I want only to repeat some basic points and to 
emphasize the differences between SU(3) and SU(2) gauge fields. The Hamiltonian for 
SU(3) gauge fields on a lattice ^ is 

H = 7^(E ^E^E^ + ^ReE[l - <^os{Ib,)]), (1) 

where g is the coupling constant, A = A/g"^, and a is the lattice spacing. In this calculation 
A = 1.1185 and a = 1. The particular choice is not important here because the final result 
can be written in a scaling fashion(lO). The electric and magnetic fields can be expressed 
in terms of the link variables Ui = exp{—i^X°'Af), which are elements of the group SU(3). 
A set of classical equations can be derived from (1): 

Ui = tY.yr^'^ui, (2) 
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= i^Et^m-t^;), (3) 



4a2 p 
where 

Vr = ^E^, (4) 

and Up is the product of all four link variables on an elementary plaquette. These equa- 
tions define a trajectory {E^{t),Ui{t)} in the classical phase space for each given initial 
condition {Ef{0), Ui{0)}. In order to study the exponential divergence of two trajectories 
we introduce the following gauge-invariant distance between two trajectories 1,2: 

Dil,2) = -^J2\^rU,p-trU2p\. (5) 

D is propotional to the absolute local difference in the magnetic energy of two different 
gauge fields. Np is the total number of plaquettes: Np = 3N^, where N is the lattice size. 
In the following calculation, we use a 10^ lattice. 

Initially we put all = to avoid the presence of local static color charges on the 
lattice. Link variables are parametrized by eight real angles a", a = 1,...,8 [|10|. Here 
there are two problems that do not exist in SU(2) simulation. First, in the case of SU(2) 
the 3 angles are used throughout the simulation. The group multiplication can be easily 
realized in this representation. But in the case of SU(3), if the 8 a's are directly used, the 
group multiplication is difficult to realize. We find that it is much faster to directly use 
the unitary 3x3 matrices in the program. Second, in SU(2) the initial magnetic energy 
put on the lattice is controlled by a single angle u. Changing the range of u changes the 
energy. But we find in our case, if we change the range of the angles in order to change 
the energy, we run the risk of restricting the initial configuration to a small subspace of 
SU(3). So here we use a different method to initialize the link variables. We start from 
a random field configuration obtained by arbitrarily selecting a's for each link. Then 



we use a heat bath algorithm [|r^ to thermalize the magnetic sub-system to temperature 
Tq. Now Tq is the only parameter to control the total energy on the lattice. A second 
configuration is chosen in the close neighborhood of the first one (-D(O) is small). This 
can be achieved by choosing slightly different a's (controlled by a parameter 6,\ 6 \<< 1) 
for each corresponding link of the second configuration with respect to that on the first 
one. 

We numerically integrate the (5) and (6). Using the measure in (8), we observe a 
similar behavior here of divergence of two trajectories as for SU(2) in Ref[l]. In Fig.l.a, 
the evolution of ln{D{t)) for Tq is shown. We see after several initial oscillations, D{t) 
increases exponentially with t and then saturates at large t due to compact nature of 
SU(3) group. We notice in the exponential increasing region the fluctuation is small for 
large Tq. This suggests the divergent property is very similar in the entire phase space 
except for some regions of small measure. We identify the slope of ln[D[t)) in the linear 
region as the largest Lyapunov exponent h. 

We have also tried a different definition of distance 

^^(^,1,2) = I JliiK? - iKf) I, (6) 
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namely the sum of the absolute value of local difference in the electric energy. The rising 
oiln[D^{t)) shown in Fig.l.b is coincident with ln{D{t)) except for the initial oscillatory 
region. This is somehow within our expectation because the chaoticity of a system is an 
intrinsic property and it shall not depend on a particular metric we choose for it. 

We have studied the wave-length dependence of the trajectory divergence. We define 
a local distance as 

L'(l,2,t,p) =1 tr[/ip-tr[/2p I, (7) 

which is a function of t and p. Here p denotes the relative position of a particular 
plaquette on the 3-dim lattice. £'(l,2,t,p) can be Fourier transformed into a function 
defined on 3-dim cubic lattice in k space which is reciprocal to the original lattice. 

D(t,k) = ^D(l,2,t,p)e^''P. (8) 
p 

To incorperate the discreteness of the lattice we define a mode spectrum as follows: 

s{t,k)^kY.\b{t,\^)\'' . (9) 

k 

which gives the relative importance of each wave length component to the divergence. In 
Fig. 2. a and Fig.2.b, 5'(0, k) and S{tf, k) for a typical run are plotted as a function of k. 
tf is the time when saturation is reached. We can see that from t — to t — tf only 
the scale of S changes. The shapes in the two cases are indistinguishable. Both of them 
are almost flat when k < 4 which indicates the randomness of the difference between the 
two configurations and decrease at larger k. So we conclude there is no evidence of any 
preference of the divergence of some specific wavelength components of the distance. The 
decreasing in k > A region is an artificial effect on a finite cubic lattice, that is, there are 
fewer lattice sites for large A; on a finite cubic lattice than we expect for an infinite one. 
Most of the big fluctuations are simply due to combinatorical properties on the discrete 
lattice, e.g. it is easier to form a 5 by the sum of squares of three integers than to form 
a 4 so we have a peak at k = \/5. Considering these artificial effects, we think the lattice 
might be too small for a conclusive analysis of this kind. 

After several runs with different initial Tq's, we get a relation between h and E shown 
in fig.3. Here E is the averaged energy on each plaquette. Using the scaling property of 
the Hamiltonian in (1), we can get a parameter invariant linear relation: 

A E 1 

ha^-g^(— )a = —g^Ea. (10) 

The number is estimated from the line in fig.3. We see the points fit nicely on the straight 
line which goes through the origin. In the classical limit, h is independent on a because 
g docs not run with a in this limit. The difference between the two points for the lowest 
To is an indication of the statistical uncertainty of the value h. The fluctuations are 
small because of the large number of degrees of freedom involved. The lack of points at 
low energy is caused by the rapid increase of time needed to thermalize a system at low 
temperature. Also as we lower the temperature, quantum effects will become more and 
more important so the correctness of a classical calculation will become doubtful. We 
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postpone this problem to the future study. To a good approximation, E is related to the 
temperature as [ p!3|] : 

16 

E=-T. (11) 

Combine (10) and (11), we get 

h = 0.54/T (12) 

It is a general result of classical non-linear dynamics that the sum of all positive 
Lyapunov exponents describes the entropy growth rate, which is roughly the inverse of 
the time scale that the system approaches thermal equilibrium (for example, see [Q). So 
the largest Lyapunov exponent leads to an estimate of the thermalization time Tg = h~^. 
If we insert the thermal coupling constant [|l| 

« = llln(.r/A)^ (l^** 

(with A = 200 MeV), we get an expression of Ts{T). As an example to show the time 
scale, at T = 400 MeV, is roughly 0.24 fm/c. So here we really see that the gluon 
system thermalizes very fast. 

Finally we want to have a short discussion on the gluon damping rate. The motiva- 
tion here is the seemingly numerical accident that the value of h is remarkably similar to 



twice the damping rate of a gluon at rest |T^, 

j(uj = 0) = 0.264/T, (14) 

in a thermal gluon system. This is a quite surprising result because these two quantities 
appear in totally different contexts and are calculated with different methods. On one 
hand, the damping rate is the imaginary part of the self energy of a quasi-particle in a 
thermal gluon system and is calculated by a sophisticated effective quantum field the- 
ory, while the Lyapunov exponent here is a classical dynamical quantity describing the 
divergent property of two classical trajectories. But though we cannot establish a direct 
relation between these two quantities in the moment, we think this similarity does not 
arise without any reason. First these two quantities, though very different from their 
contexts, both describe how fast a non-equilibrated gluon system approaches thermal 
equilibrium. The relevance of h is clear from the above discussion. The connection of 
7(0;) appears clearly [|1^] in 

/K^) = ^IJ7^ + c(a^)e-^^^^^ (15) 

where / is time dependent distribution function, the first term on the l.h.s is the equili- 
brated distribution. We see 7 appears in the decaying term. Second, we want to prove 
that the gluon damping rate is basically a quantity of semi-classical origin. Since 7(0;) 
is a smooth function, it is sufficient to prove the statement for a gluon at rest for which 



case explicit quantum field calculation is given in Ref.|lJ|. So in the following, we will 
concentrate only on the gluon at rest. 
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In Ref.[n|, the damping rate is obtained from the soft part of the thermal self energy 
loops of a gluon at rest (Fig. 4. a). Equivalently, it can also be obtained by calculating cross 
sections from the diagrams resulting by cutting the self energy loops, multiplying by the 
thermal distribution functions and integrating over phase space [|l5l. Or physically, we 
can first calculate the decay rate 7^ and the production rate 7^ of a gluon at rest in a 
thermal gluon system and then substract them to get 7 = 7d — 7p- From the cutting of the 
self energy loop, we can see that the contributions to 7 come from Rutherford scattering 
(Fig.4.b) and bremsstrahlung processes (Fig.4.c). Both of them need only to be calculated 
on the tree level to get the leading order contributions to 7 (of the order of g^T) , provided 
the effective Green's functions are used for the soft vertices and soft propagators appearing 



in the diagrams. This latter requirement is equivalent to resummation [|16 



On the other hand, we know generally the result of a tree diagram in a massless 



field theory is of classical limit [|T3- Now in our finite temperature case quantum effects 
can come from two directions. First the effective vertices and propagators might be of 
quantum origin. Secondly we must use the quantum thermal distribution function in 
order to get a finite result of thermal contribution. 

First we consider the effective Green's functions. We want to prove that these 
functions are of semi-classical origin in the high temperature limit. It is shown ||T^ that 
in the high temperature limit the leading order modification to a soft propogator or 
vertex is obtained by inserting a hard thermal loop (HTL) in the bare quantity (Fig.S.a 
and Fig.S.c). The hard lines themselves need not be modified if we only require leading 
order solutions. The thermal contributions correspond to diagrams with a cut through 
one internal line of the hard thermal loop at any place (put the corresponding line on 
mass shell). This can be seen in the real time formalism of finite temperature field 
theory. In this formalism, a thermal propagator is the sum of the vacuum propagator 
and another term due to finite temperature D ~ n{uj)5{p'^ — rn^), where n{uj) is the 
thermal distribution function and the 5-function puts the particle on mass shell. The 
multiplication of the vacuum part of the propagators around the loop just gives the 
vacuum contribution which leads to renormalization. The thermal contribution comes 
from terms each of which has a (5-function from the second term of a propagator. The 
cutting procedure changes the original loop diagram into a tree level diagram (Fig.S.b and 
Fig.S.d) with thermal distribution functions attached to the two new external legs. So we 
see that the modification to a soft propagator or a vertex due to a hard thermal loop is 
at tree level, i.e., it is of semi-classical origin. The quantum effects come only from the 
Bose distribution functions. To suggest a possible method to actually obtain these HTL's 
classically, here we mention that the polarization effect (2-point HTL) in a QED plasma 
was obtained long ago within classical kinetic theory [0 . Recently the kinetic formalism 
has been used to calculate all hard thermal loops ||T9|| . 

Now in principle we can start from these effective quantities and construct a classical 
perturbative method. Then we can use it to calculate bremsstrahlung processes and 
Rutherford scattering and finally obtain the gluon damping rate to order g^T. 

We have proved that 7(0) is a semi-classical quantity. We now further show that it 
is really classical. We have seen that the quantum effects come into 7(0) only through 
a trivial way, namely, via the Bose distribution function. If we want to use the classical 
phase space distribution n = T/tu, which corresponds to the high temperature limit of 
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Bose distribution, in the calculation, we must use a ultraviolet cutoff (lattice size a) in 
order to avoid divergences. Then all the HTL's may depend on a. Now to see how a comes 
into a quantity, we shall express this quantity in terms of classical variables: g"^ = Anas/h, 
T and a. Here Qc is so defined that h does not appear explicitly in classical Yang-Mills 
equations. For example, the plasma frequency will be ujp = Qc^jT/a, which diverges when 
a goes to zero. Though Up can be calculated semi-classically, it is not well defined at 
strict classical limit. On the other hand 7(0) ~ giT is independent on this cutoff a in the 
leading order. This proves that to the leading order 7(0) is a purely classical quantity. 

To conclude, we have shown that classical SU(3) gauge theory is chaotic. This 
provides a classical explanation for why the gluon system thermalizes rapidly. 
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Figure Captions 



Fig.l Divergence of two classical trajectories in two different measures D{t) (a) and D^(t) 
(b). 

Fig. 2 Wave-length spectrum of the local distance t — 0(a) and t — tf{\)). 
Fig. 3 Lyapunov exponent as a function of energy. 

Fig.4 (a) Soft part of the self energy loop for a gluon at rest, (b) Rutherford scattering 
and (c) bremsstrahlung process resulting from the cutting of (a). The solid lines 
are for hard gluons with momenta of the order of T and the dotted lines are for soft 
ones of the order gT. The blobs denote effective Green's functions resulting from 
inserting hard thermal loops in diagrams with soft external gluon lines. 

Fig. 5 (a) HTL contribution to the self energy of a soft gluon, (b) tree level diagram from 
cutting of (a), (c) HTL contribution to the soft three gluon vertex, (d) tree level 
diagram from cutting of (c). Lines have same meaning as that in Fig.4. 
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